import numpy as np

# 参数设置
U = 1.0       # 来流速度
a = 1.0       # 圆柱半径

# 理论推导部分：
# 当圆柱绕流的环量 Γ 超过临界值 4πUa 时，驻点会脱离圆柱表面。
# 流函数ψ表达式：
# ψ = U (r - a²/r) sinθ + (Γ/(4π)) ln r
# 径向速度 u_r 和切向速度 u_θ：
# u_r = U (1 - a²/r²) cosθ
# u_θ = -U (1 + a²/r²) sinθ + Γ/(2πr)
# 驻点条件：u_r = 0 和 u_θ = 0
# 在圆柱表面 r = a 处，驻点条件变为：
# -U (1 + 1) sinθ_s + Γ/(2πa) = 0
# 解得：sinθ_s = Γ/(4πUa)
# 当 Γ = 4πUa 时，sinθ_s = 1，即 θ_s = π/2 或 3π/2。
# 当 Γ > 4πUa 时，sinθ_s > 1，无解，说明驻点脱离圆柱表面。

# 计算临界环量
Gamma_critical = 4 * np.pi * U * a
print(f"临界环量 Γ_critical = {Gamma_critical:.2f}")

# 工程修正：引入安全系数
n_min = 1.5
n_max = 2.0
Gamma_safe_min = Gamma_critical / n_max
Gamma_safe_max = Gamma_critical / n_min
print(f"安全环量范围 Γ_safe ∈ [{Gamma_safe_min:.2f}, {Gamma_safe_max:.2f}]")

# 参数扫描范围
Gamma_min = 0.0
Gamma_max = 6 * np.pi * U * a
print(f"参数扫描范围 Γ ∈ [{Gamma_min:.2f}, {Gamma_max:.2f}]")

# 输出安全环量范围
print("\n安全环量范围确定结果：")
print(f"当 Γ > {Gamma_critical:.2f} 时，驻点脱离圆柱表面，流动可能分离。")
print(f"建议的安全环量范围为 [{Gamma_safe_min:.2f}, {Gamma_safe_max:.2f}]，其中 n ∈ [1.5, 2.0]。")